#Binf739-2: Gene Regulatory network Final Project
#Team Member:
#Nuttachat Wisittipanit
#Saeed Khoshnevis

import string
from numpy import *
from scipy import *
import copy


#######################
#####Main Program######
#######################

#Get the network using BFS (Breadth-First Search Algorithm)
def getnetworkbreath(fa,fx,target='CG9885',startW=0,endW=30,mcase=9,proc=1.94,kn=.8,loopdepth=2,seed=['z','abd-A','dl','twi','brk','bap','ctBp','nej','pnr','tin'],outall='outall.txt',outseed='outseed.txt',outmseed='outmseed.txt'):
    es=run(fa)
    gg=runx(fx,target)
    bg=gobt(es,target,startW,endW,kn,loopdepth,mcase,proc)
    nre,ncmoreseed,ncseed,bign,nwseed,ngseed=clpath(gg,es,bg,seed)
    adjall,nodeall=adjm(bg,bign)
    adjseed,nodeseed=adjmx(nre,ncseed)
    adjmoreseed,ncmoreseed=adjmx(nre,ncmoreseed)
    Save(outall,adjall)
    Save(outmseed,adjmoreseed)
    Save(outseed,adjseed)
    return nodeall,adjall,ncmoreseed,adjmoreseed,nodeseed,adjseed,bg,nwseed,ngseed

#Get the network using DFS (Depth-First Search Algorithm)
def getnetworkdepth(fa,fx,target='CG9885',startW=0,endW=30,mcase=9,proc=1.94,kn=.8,loopdepth=2,seed=['z','abd-A','dl','twi','brk','bap','ctBp','nej','pnr','tin'],outall='outall.txt',outseed='outseed.txt',outmseed='outmseed.txt'):
    es=run(fa)
    gg=runx(fx,target)
    bg=depth(es,target,startW,endW,kn,loopdepth,mcase,proc)
    nre,ncmoreseed,ncseed,bign,nwseed,ngseed=clpath(gg,es,bg,seed)
    adjall,nodeall=adjmx(bg,bign)
    adjseed,nodeseed=adjmx(nre,ncseed)
    adjmoreseed,ncmoreseed=adjmx(nre,ncmoreseed)
    Save(outall,adjall)
    Save(outmseed,adjmoreseed)
    Save(outseed,adjseed)
    return nodeall,adjall,ncmoreseed,adjmoreseed,nodeseed,adjseed,bg,nwseed,ngseed

#######################
#######################


#######################
#####Get Dataset#######
#######################

def teaa(a=3,b=2,c=5):
    print a,b,c
    
def readfile(fn):
    fp = open(fn)
    data = fp.read()
    fp.close()
    return data

def Mlist(data):
    n   = data.count('\n')
    c   = 1
    l   = []
    for i in range (n+1):
        iN  = data.find('\n',c)
        #if i==9:
        #    print c
        if i != n:
            l.append(data[c-1:iN])
        else:
            l.append(data[c-1:len(data)])
        c   = iN+1
        #print i,iN,c
    #print len(l)
    return l

def clean(l1):
    e1  = []
    for i in range (len(l1)):
        e1.append(l1[i].split())
    e1.pop()
    return e1

def cs(l1):
    e1=[]
    for i in range (len(l1)):
        e1.append(l1[i].split())
    #e1=array(e1)
    for i in range(len(e1)):
	for j in range(1,len(e1[i])):
                #print j
		e1[i][j]=e1[i][j]
    e1.pop()
    return e1

##############################
#####Data Normalization#######
##############################

def toStandard(ListP):
    Slist=ListP[:]
    n=len(ListP)
    for i in range (n):
        temp=zeros(31,float)
        for j in range(31):
            temp[j]=eval(ListP[i][j+1])
        sd=0.0
        sd=temp.std()
        m=0.0
        m=temp.mean()
        for k in range(31):
            temp1=(eval(ListP[i][k+1])-m)/sd
            temp1=round(temp1)
            if temp1==0:
                Slist[i][k+1]=0.0
            elif temp1 <0:
                Slist[i][k+1]=-1.0
            elif temp1 >0:
                Slist[i][k+1]=1.0
    return Slist
def run(fn):
    data    = readfile(fn)
    l1      = Mlist(data)
    e1      = clean(l1)
    es1     = toStandard(e1)
    return es1
def runt(fn):
    data = readfile(fn)
    l1=Mlist(data)
    e1=cs(l1)
    return e1
def runx(fx,target):
    data=readfile(fx)
    m=Mlist(data)
    n=cs(m)
    nln=[]
    for i in range(len(n)):
        for j in range(len(n[i])):
            if (n[i][j]!=target) and (n[i][j] not in nln):
                nln.append(n[i][j])
    #print len(n)
    return nln

##################################
###########Get Pathways###########
##################################

#ITFF (Immediate Factor Finding Algorithm)
def getimd(md,tgdummy,nodedummy,startW,endW,mcase,kn,proc):
    #print kn
    mdcc=copy.deepcopy(md)
    whl=len(mdcc)
    tgdcc=copy.deepcopy(tgdummy)
    ndcc=copy.deepcopy(nodedummy)
    subnode=[]
    
    eprob=[]
    nprob=[]
    cprob=[]
    N=len(md)
    n=0
    c=0
    #print tgdcc,len(mdcc),ndcc,startW,endW,kn,mcase,proc,
    for i in range(whl):
        
        if mdcc[i][0] not in ndcc:
            #if i==0:
            #    print md[i]
            #print mdcc[i][0],
            temp1=getDependantProb(tgdcc,mdcc[i],startW, endW)
            if temp1 > 0:
                temp=getEnhanceProb(tgdcc,mdcc[i],i,startW, endW,mcase)
                #print kn
                if temp[0]> kn :
                    eprob.append(temp)
                    #print len(eprob)
                temp2=getCounterProb(tgdcc,mdcc[i],i,startW, endW,mcase)
                #print temp2,
                if temp2[0]> kn :
                    cprob.append(temp2)
            else:
                c+=1
    #print len(eprob),len(cprob)           
    total=sumProb(eprob,cprob)
    #print total,
    child=getChildren(total,proc)
    ch=getname(md,child)
    #print len(total),child,ch,
    for i in range(len(ch)):
        node=[tgdcc[0],ch[i]]
        nodetemp=copy.deepcopy(node)
        subnode.append(nodetemp)
        node=[]
    
    return subnode

def getds(nm,es):
    for l in range(len(es)):
        if nm==es[l][0]:
            ds=es[l]
    return ds   

#Main funtion for BFS
def gobt(md,target,startW,endW,mcase,loopdepth,kn,proc):
    #print kn
    node=[]
    bignode=[]
    nodedummy=[]
    node.append(target)
    bignode.append([target])
    nodedummy.append(target)

    nodes=copy.deepcopy(node)
    chknode=len(nodes)
    loopchk=0
    rc=1
    
    while rc:
        node=[]
        for i in range(chknode):
            #print i,nodes[i]
            tgdummy=getds(nodes[i],md)
            #print tgdummy[0]
            subnode=getimd(md,tgdummy,nodedummy,startW,endW,kn,mcase,proc)
            #print subnode
            if len(subnode)==0:
                continue
            #print tgdummy,subnode
            for j in range(len(subnode)):
                if subnode[j][-1] not in node:
                    node.append(subnode[j][-1])            
            bignode=conn(bignode,subnode)
        
        nodes=copy.deepcopy(node)
        chknode=len(nodes)
        
        for i in range(chknode):
            nodedummy.append(nodes[i])

        loopchk+=1
        if loopchk==loopdepth or chknode==0:
            rc=0
    
    return bignode

#Main function for DFS
def depth(md,target,startW,endW,kn,loopdepth,mcase,proc):
 
    node=[]
    bignode=[]
    histnode=[]
    
    tg=getds(target,md)
    #print tg
    #sd=getds(seed,md)

    node.append(target)
    nodedummy=node[:]
    nodes=node[:]
    histnode.append(tg)

    loopchk=0
    chknode=len(node)
    
    rc=1
    oj=1
    ckmd=len(md)
    ckmdp=ckmd

    hnd=histnode[:]
    tgdummy=hnd.pop()
    chktg=[]
    
    ct=0
    #print nodes
    while rc:
        while oj:
            #print tgdummy
            if (md[ckmd-ckmdp][0] not in nodedummy):
                #print md[ckmd-ckmdp][0]
                temp1=getDependantProb(tgdummy,md[ckmd-ckmdp],startW, endW)
                if temp1>0:
                    temp=getEnhanceProb(tgdummy,md[ckmd-ckmdp],(ckmd-ckmdp),startW, endW,mcase)
                    temp2=getCounterProb(tgdummy,md[ckmd-ckmdp],(ckmd-ckmdp),startW, endW,mcase)
                    #print temp,temp2
                #print kn,loopdepth,proc,mcase
                if  (temp[0]> kn) and (temp2[0]> kn) and ((temp[0]+temp2[0])> proc):
                    #print 'kkk',
                    node.append(md[ckmd-ckmdp][0])
                    nodedummy.append(md[ckmd-ckmdp][0])
                    nodes=node[:]
                    if tgdummy[0] not in chktg:
                        chktg.append(tgdummy[0])
                    #print chktg
            if (chknode + 1) == len(node):
                chknode = len(node)
                histnode.append(md[ckmd-ckmdp])
                hnd     = histnode[:]
                tgdummy = hnd.pop()
                loopchk += 1
                ckmdp   = ckmd
                
            else:
                ckmdp-=1
            if ckmdp==0 or loopchk==loopdepth:
                oj=0
                
        if ckmdp!=0:
            bignode.append(nodes)
            ckmdp=ckmd
            histnode.pop()
            hnd=histnode[:]
            try:
                tgdummy=hnd.pop()
            except IndexError:
                tgdummy=[]
            oj=1
            loopchk-=1
            #print 'here',node
            node.pop()
            chknode=len(node)
            #print bignode
        else:
            if tgdummy[0] not in chktg:
                bignode.append(nodes)
            oj=1
            ckmdp=ckmd
            loopchk-=1
            histnode.pop()
            hnd=histnode[:]
            try:
                tgdummy=hnd.pop()
            except IndexError:
                tgdummy=[]
            #print node
            node.pop()
            chknode=len(node)
        ct+=1
        print ct,
        if chknode==0:
            rc=0
    
    return bignode

def conn(bignode,subnode):
    fst=copy.deepcopy(bignode)
    alt=copy.deepcopy(subnode)
    upl=len(fst)
    lwl=len(alt)
    fstx=[]
    for i in range(upl):
        if fst[i][-1]==alt[0][0]:
            for j in range(lwl):
                buff=copy.deepcopy(fst[i])
                buff.append(alt[j][1])
                fstx.append(buff)
        else:
            fstx.append(fst[i])
    return fstx

def getChildren(total,proc):
    child=[]
    for i in range(len(total)):
        #if float(total[i][1]) >=max(array(total)[:,1]):
        if float(total[i][1]) > proc:
	    child.append(total[i])
    return child
def getname(md,child):
    nm=[]
    for i in range(len(child)):
        nm.append(md[child[i][0]][0])
    return nm
def getEnhanceProb(tg,mdi,index,startW,endW,mcase):
    cnt=cnts=0
    #print startW,endW,
    for i in range(startW,endW):
        if  (mdi[i+1]> 0 and  tg[i+2]> 0):
            cnt +=1
        if  (mdi[i+1]>0):
            cnts +=1
    if cnts==0 or cnts < mcase:
        prob=0
    else:
        prob=float(cnt)/(cnts)
    #print cnt,cnts,
    return [prob,index]

def getDependantProb(tg,mdi,startW, endW):
    cnt=cnts=dnts=0
    #print tg,len(mdi),startW,endW
    for i in range(startW, endW):
        if  (mdi[i+1]> 0 and  tg[i+2]> 0):
            cnt +=1
        if  (mdi[i+1]>0):
            cnts +=1
        if  (tg[i+2]>0):
            dnts +=1
    if cnts==0 :
        dprob=0
    elif abs((float(cnt)/(cnts))-(float(dnts)/(endW-startW))) > 0:
        dprob=1
    else:
        dprob=0
    #print cnt,cnts,dnts,
    #print dprob,
    return dprob

def getCounterProb(tg,mdi,index,startW,endW,mcase):
    cnt=cnts=0
    for i in range(startW,endW):
        if  (mdi[i+1]<= 0 and  tg[i+2]<= 0):
            cnt +=1
        if  (mdi[i+1] <= 0):
            cnts +=1
    if cnts==0 or cnts < mcase:
        prob=0
    else:
        prob=float(cnt)/(cnts)
    #print cnt,cnts,
    return [prob,index]

def sumProb(eprob,cprob):
    total=[]
    for i in range(len(eprob)):
        temp=eprob[i][1]
        for j in range(len(cprob)):
            if temp==cprob[j][1]:
                total.append([cprob[j][1],(cprob[j][0]+eprob[i][0])])
                break
    return total

#Get nodes
def clpath(gg,es,bg,seed):
    esg=copy.deepcopy(es)
    cseed=copy.deepcopy(seed)
    nbg = copy.deepcopy( bg )
    ggg=copy.deepcopy(gg)
    newbg=[]
    ln=-1
    
    nlseed=[]
    nwseed=[]
    ngseed=[]
    
    for i in range(len(esg)):
        if (esg[i][0] in ggg) and (esg[i][0] not in nlseed):
            nlseed.append(esg[i][0])        
    for i in range(len(nbg)):
        for j in range(len(nbg[i])):
            if (nbg[i][j] in nlseed) and (nbg[i][j] not in nwseed):
                nwseed.append(nbg[i][j])
    for i in range(len(nlseed)):
        if (nlseed[i] not in nwseed) and (nlseed[i] not in ngseed):
            ngseed.append(nlseed[i])
    

    for i in range(len(nbg)):
        for j in range(len(nbg[i])):
            if nbg[i][j] in cseed:
                newbg.append(nbg[i])
                break
    nng=copy.deepcopy(newbg)
    for i in range(len(nng)):
        ok=1
        while ok:
            if nng[i][ln] not in cseed:
                nng[i].pop()
            else:
                ok=0
    nna=copy.deepcopy(nng)
    #print nna
    nresult=[]
    #nresult.append(nna[0])
    for i in range(len(nna)):
        if nna[i] not in nresult:
            nresult.append(nna[i])
    ncr=[]
    ncr.append(nbg[0][0])
    for i in range(len(nresult)):
        for j in range(len(nresult[i])):
            if nresult[i][j] not in ncr:
                ncr.append(nresult[i][j])
    #print len(nbg)
    ncrr=[]
    ncrr.append(nbg[0][0])
    for i in range(len(ncr)):
        if (ncr[i] in cseed) and (ncr[i] not in ncrr):
            ncrr.append(ncr[i])
    
    bignode=[]
    for i in range(len(nbg)):
        for j in range(len(nbg[i])):
            if nbg[i][j] not in bignode:
                bignode.append(nbg[i][j])
     
    return nresult,ncr,ncrr,bignode,nwseed,ngseed



###################################
######Create Adjacency Matrix######
###################################

def adjm(pathway,node):
    print 'hereadj',
    d=len(node)
    nodex=copy.deepcopy(node)
    pw=copy.deepcopy(pathway)
    adm=zeros((d,d),int)
    q=0
    for i in range(d):
        for j in range(i):
            a=nodex[i]
            b=nodex[j]
            for k in range(len(pw)):
                for l in range(len(pw[k])-1):
                    if (pw[k][l]==a and pw[k][l+1]==b) or (pw[k][l]==b and pw[k][l+1]==a):
                        adm[i][j]=1
                        adm[j][i]=1
                        q=1
                        break
                if q==1:
                    q=0
                    break
                    
                
    return adm,nodex
                
def adjmx(pathway,node):
    print 'hereadjx',
    d=len(node)
    nodex=copy.deepcopy(node)
    pw=copy.deepcopy(pathway)
    adm=zeros((d,d),int)
    q=0
    for i in range(d):
        #print '1st'
        for j in range(i):
            #print '2nd'
            a=nodex[i]
            b=nodex[j]
            for k in range(len(pw)):
                try:
                    ax=pw[k].index(a)
                    bx=pw[k].index(b)
                    #print ax,bx
                    if (ax == bx+1) or (ax+1==bx):
                        adm[i][j]=1
                        adm[j][i]=1
                        print 'breakfux'
                        break
                except ValueError:
                    #print 'yeh'
                    continue

    return adm,nodex

#####################################
######Export the Adjacency File######
#####################################
rf='c:/program files/r/r-2.4.1/'
def Save( name, data ):
# the number of rows
    if len(data.shape)>1:
        (nc,nr) = data.shape
    else :
        nc = data.shape[0]
        nr = 1
    fp = open( rf+name, 'w')
    if nr > 1 :
        for i in range(0,nc):
            for j in range(0,nr):
                #print j,
                fp.write( str( data[i,j] ))
                fp.write(' ')
            fp.write(' ')
            fp.write('\n')
        fp.write('\n')
    else:
        for i in range(0,nc):
            fp.write( str(data[i]) )
            fp.write(' ')
        fp.write( '\n' )
    fp.close()

    



